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Abstract. We performed Gibbs Ensemble Molecular Dynamics simulations to 
determine the Liquid- vapor coexistence curves of metals modeled with a modified form 
of inter-particle pair potential. The parameters of the potential are obtained by fitting 
the cold curve (Energy/atom Vs volume curve at OK) obtained from the potential to 
that obtained from ab-initio calculations. Simulations are done for Aluminum, Copper, 
Sodium and Potassium and the results are analyzed. We find that the present results 
improve significantly over those obtained from Morse potentials (J.K. Singh et. al., 
Fuid Phase Equilibria 248(2006)). 



1. Introduction 

Equation of State(EOS) of materials for wide ranges of temperatures, pressures and 
densities is an essential requirement of hydrodynamic simulations which find applications 
in astrophysics, high energy density physics etc. First principles based density functional 
theory(DFT) calculations provide accurate description of thermodynamic properties of 
materials at zero temperature. However, to do a non-zero temperature calculation, 
particularly in the fluid regime, DFT has to be coupled with either Molecular Dynamics 
(MD) or Monte-Carlo (MC) methods. This enormously increases the computational cost 
and the method becomes impractical for large scale computations which are required 
for many practical purposes. Thus theoretical modeling of EOS of materials in the 
expansion phase based on classical statistical mechanical methods still plays a major 
role in generating thermodynamic data of materials. Basic input for theoretical modeling 
of EOS is the inter-particle interaction. 

Presently, various forms of embedded atom model(EAM)[U [2] potentials are 
available in the literature. However, a recent simulation study of liquid-vapor phase 
diagram using a particular form of EAM showed that the potential even though is 
good in predicting solid state properties, is not accurate enough in predicting the liquid 
vapor phase diagram [3]. Besides that, all the presently available theoretical models 
for EOS are restricted to pairwise inter-particle potentials. A way to obtain effective 
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pair potential from EAM potential has been developed by Foiles[l] based on Taylor 
series expansion of the embedding term around equilibrium solid density. However, the 
method is found to be working only for densities close to melting. Thus obtaining pair 
potentials which can predict accurately their thermodynamic properties for wide ranges 
of temperatures and pressures is an important requirement. 

In the present paper, we emphasize on metals. Two model pair potentials have been 
widely used to describe the inter-particle interactions in metals. They are the generalized 
Lennard- Jones (GLJ) and the Morse potentials. The parameters of Morse potential for 
various cubic metals have been obtained by Lincoln et. al. [5] by using the Lattice 
parameter, bulk modulus and cohesive energy data. Lincoln et. al. have accounted for 
interaction between particles up to ninth neighbor shell in obtaining the parameters. 
J.K. Singh et. al. [6] have obtained the liquid-vapor phase diagrams for various metals 
using the Morse potential as parametrized by Lincoln et. al. from classical MD. Their 
results show that there is enormous deviation between the simulation and experimental 
values in the case of alkali metals although the results show a reasonable agreement 
with literature data in the case of Aluminum , Copper and Gold. On the other hand, 
there is a recent parametrization of the GLJ potential with some modification by Sun 
Jiuxun[7]. The author showed that Pressure Vs volume curve at OK generated from 
this potential by taking only nearest neighbor interaction is in good agreement with 
experimental results. However, we observed that the energy Vs volume curve is not 
matching with the ab-initio data and the cohesive energy is highly erroneous in some 
cases jH]- Apart from this, the excess internal energy per particle calculated from the 
GLJ potential using parameters obtained by Sun Jiuxin is found to be diverging for most 
of the materials and classical molecular dynamics was not possible with that potential. 

These problems might have araised because of attempting to fit the potential with 
only two or three experimentally measured physical quantities like cohesive energy, 
bulk modulus and its pressure derivative, which are just equilibrium properties. The 
accuracy of the cold curve away from the equilibrium particularly in the expansion 
region is unknown . However, accuracy of the cold curve in the expansion region might 
play a role in determining the accuracy of the calculations of thermodynamic properties 
of corresponding fluid. On the other hand, accurate and fast DFT calculations of the 
cold curve are possible with the present day available codes. Thus we use the cold curve 
obtained from ab-initio calculations to obtain parameters of the potential. We fit the 
cold curve obtained from a given form of potential to that obtained from ab-initio data. 
This gives a clear idea of the accuracy of the parametrized cold curve form far from the 
equilibrium. When this procedure is done with Morse and GLJ potentials, we observed 
that the cold curves obtained from them are deviating from ab-initio data away from 
the equilibrium particularly in the case of alkali metals. Hence we propose a modified 
form of empirical potential which we expect to be accurate enough even far from the 
equilibrium for a wide class of metals. 

Thus in the present paper, we describe an alternative way of obtaining parameters 
of empirical potentials and also a modified form of empirical potential. To test the 
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accuracy of the potentials so obtained, we obtained liquid vapor coexistence points 
for various temperatures for Aluminum, Copper, Sodium and Potassium using Gibbs 
ensemble molecular dynamics(GEMD) with the modified potential. In Section II, we 
describe the way we obtained parameters of the potentials using ab-initio data and the 
cold curves (Envergy Vs volume curves) from different potentials for various metals are 
compared and contrasted. In sec III, the GEMD method we adopted is explained in 
brief and liquid-vapor phase diagrams we obtained for Aluminum, Copper, Sodium and 
Potassium are shown and results are analyzed. The paper is concluded in Section IV. 



2. Pair Potential Models for Metals 



The GLJ potential is given by 
u L j(r) 



mi 



. , "-I - (1) 
mi — rii i v r J \r J . 

Jiuxun[7] proposed a relation between m\ and n\ so that only one is independent. The 
relations are m\ = 6n — 3 and n\ = 3n — 3. Taking his modification into account, the 
potential becomes, 

ulj{t) = — (3n - 3) — - (6n - 3) — 
on |_ \r / \ r J 

In the present work we consider parameter. The Morse potential is 

UM (r) = e ( e -Mr-ro) _ 2e -a(r-ro)) 

We propose a new form of pair potential given by 

-2«(r/ro-l) ( r A\ W _ 2e -a(r/ro-l) ' 

\r J \ r 

Attractive part of the above potential is chosen inspired by the screening of ions by 
electrons in metals. Repulsive part is chosen ad-hoc as per mathematical convenience. 
It can be seen that by putting = in above potential we recover the Morse potential 
and by putting a = 0, it becomes a Lennard- Jones type potential. Since the electron- 
ion screening term is represented by an exponential term multiplied by the coulomb 
term (at-least in the semi-classical picture), we expect that a similar form with the 
exponents of the exponential and the r in the potential being adjusted, would represent 
the interactions in a better way. In the case of GLJ potential, ro,n and e are the 
parameters. For the Morse potential, e, a and r are parameters and for the modified 
potential, e, a, (3, r are the parameters. For all the potentials, the potential is minimum 
at r and the well depth is e. 



u s r 



(2) 



(3) 



(4) 



2.1. Obtaining the Potential Parameters 

We used Energy Vs volume data obtained from ab-initio calculations to obtain the 
parameters of the potentials. Ab-initio calculations were done using VASP[9] software. 
PAW potentials and Monkhorst-Pack grid in reciprocal space have been used. The 
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number of points in reciprocal space have been adjusted inspecting the convergence of 
the total energy for each case. Also, the energy of free atom has been calculated by 
taking a large volume per atom and has been subtracted from the total energy per 
particle so that cohesive energy can be calculated directly from the curve. The cold 
curve i.e., the zero Kelvin isotherm for each model potential is written as follows. 

m £ 

U = Ey^( fl i) ( 5 ) 

i=l 

where x can be LJ, M or s. U is the energy per particle, is the distance of the i th 
neighbor from a particle situated at origin. Si is the number of i th neighbors. In the 
present work we have accounted interaction up to (m =)10 th neighbor. Now cij is related 
to volume per atom as follows: 

In the case of a fee solid, Oj = yia% and a% — a/j. Where a is the lattice parameter 
and 7 is the structural constant. For fee solids, volume per atom V = a 3 /4 and 7 is 
equal to a/2. For simplicity, we write r as (4Vo) 1 ^ 3 /7 without loss of generality. Thus 
using this information, Eq.([5]) can be written in terms of volume per atom V. The 
equation then becomes 

m r 

U = J2fu x (V/V ) (6) 
1=1 

For example, using u s (r) in U, 
m Se 

i=i z 

In the case of a bcc solid, 7 = 2/\/3 and V = a 3 /2. However, the a« do not hold a 
general relation with a as for fee solids and have to be carefully calculated. In this case 
r is chosen as (2V ) 1//3 . Writing <Zj = c^a where di is to be calculated using the crystal 
structure, equation corresponding to u s {f) is 

m Se 

h 2 

The parameters of the potential are obtained by fitting the cold curve Eq.(j6]) to the 
ab-initio cold curve. 

The fitted zero Kelvin isotherms obtained from various potentials for Aluminum, 
Copper, Sodium and Potassium are shown in Figs.([T]|3|). From the figures it can be 
seen that the cold curve from GLJ potential with Jiuxun's modification is quiet off from 
the ab-initio data away from the equilibrium in all the cases and in both compression 
and expansion phases. The cold curve obtained from Morse potential with parameters 
given by Lincoln et. al. (which we denote by UMLiniV)), in the case of Aluminum has 
slight deviation from ab-initio data near the equilibrium and away from equilibrium also 
as the cohesive energy predicted by ab-initio calculations and that used by Lincoln et. 
al. differed slightly. On the other hand, the cold curve from Morse potential fitted to 
the ab-initio data (denoted as Um(V)) and that from modified potential (U S (V)) had 



(7) 



-2a(di(y/y ) 1/3 -i) 



d. 



2/3 
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Figure 1. Zero Kelvin isotherm for Aluminum obtained using various potentials by 
fitting to ab-initio data. (Dotted Line: cold curve from GLJ potential), (Double dots: 
Cold curve from Morse potential parametrized by Lincoln et. al.[5]), (Dashes: Cold 
curve from Morse potential fitted to ab-initio data), (Solid line: Cold curve from 
modified potential), (circles: ab-initio data). Panel (a) is for expansion phase and 
panel (b) is for compression phase. 




Figure 2. Zero Kelvin isotherm for Copper. Description same as in Fig.([T|) 
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Figure 3. Zero Kelvin isotherm for Sodium. Description same as in Fig.([T]) 




Figure 4. Zero Kelvin isotherm for Potassium. Description same as in Fig.([T]) 
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Figure 5. Modified pair potentials for Aluminum, Copper, Sodium and Potassium. 



same accuracy except for a slight deviation of Um{V) from the ab-initio data away from 
equilibrium in the expansion phase. In the case of Copper, Umuu(V), Um(V) and U S (V) 
are equally accurate in the expansion phase. But in the compression phase, Umiau^Y) 
had a significant deviation from ab-initio data. For sodium and potassium, cold curves 
from ab-intio calculations and that given by UMiAnty) are deviating significantly This 
could be because of improved accuracy of the present ab-initio calculations in predicting 
the cohesive energies. Also, the Um{V) deviated significantly from the ab-initio data in 
the expansion region far from the equilibrium for both sodium and potassium. Whereas 
in the compression region, Um{V) is matching reasonably well with ab-initio data for 
Sodium and is deviating from ab-initio data for Potassium. On the other hand, U S {V) 
matched with reasonably well with ab-initio cold curve in both the cases. 

From the above discussion, U S (V) obtained from u s (r) matched with the ab-initio 
cold curve better than those obtained with other potentials for all the cases shown far 
away from equilibrium, particularly in the expansion region. Of course, it has to be 
noted Um(V) is accurate enough for Aluminum and Copper. Thus we expect that the 
potential u s (r) may give an improved description of the metals in the liquid region. 
u 8 (r) for Aluminum, Copper, Sodium and potassium is shown in Fig.Q Parameters for 
the potential u s (r) obtained using the procedure described above are listed in Table. 1. 
It can be seen that for Sodium and Potassium, the parameter /3 becomes negative which 
will make the potential turn down and go to zero close to origin which is unphysical. 
Thus it has to be cutoff at an appropriate point close to the origin. We found that 0.25r 
and 0.4r can be the cut-off points for Sodium and Potassium respectively. It is assumed 
that for distances smaller than these, the potential is kept constant and equal to that at 
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Table 1. Parameters for u s (r) 



metal 


e(eV) 


r (A°) 


a 


P 


Al 


0.220 


3.568 


2.499 


0.7808 


Cu 


0.321 


2.881 


3.095 


0.792 


Na 


0.086 


4.567 


3.968 


-0.5573 


K 


0.088 


5.425 


5.172 


-1.439 



the cutoff point. Clearly, this would not effect the results as the probability of finding 
an atom inside the repulsive core is negligible. In order to test the accuracy of modified 
potential u s (r) in the fluid phase, we obtained liquid- vapor phase diagrams(LVPDs) for 
all these metals using Gibbs ensemble molecular dynamics. 

3. Liquid- Vapor Phase Equilibria 

We use GEMD([T0]-[15]) method to obtain the liquid- Vapor phase coexistence curves. 
The basic idea of GEMD is to simulate the conditions of liquid- vapor coexistence. 
The system contains two simulation boxes. Total number of particles in the system 
is kept constant. However, exchange of particles between boxes is allowed. One box 
is assumed to be situated in an infinite medium of homogeneous liquid at a (given) 
constant temperature and the other is assumed to be situated in an infinite medium of 
homogeneous vapor at the same temperature. Since the idea is to get the thermodynamic 
properties of a macroscopic system, the interface effects are neglected. The coexistence 
conditions are simulated by evolving the boxes in such a way that they have same 
temperature, pressure and chemical potential after equilibration as required by the Gibbs 
phase rule. Initially, each box is given a guess density (i.e., no. of particles and volume). 
Equilibration would be faster if the guess densities are closer to the coexisting liquid 
and vapor densities. Periodic boundary conditions are applied to each box to ensure 
that they represent the bulk coexisting phases. Temperature fluctuations in each box 
are controlled by a Berendsen thermostat [16J so that they reach the given temperature. 
Pressures in both the boxes are equalized by controlling the volume fluctuations using 
a Berendsen barostat. This is done by adjusting the volume of each box such that the 
instantaneous pressure in one box becomes equal to the instantaneous pressure in the 
other. However the total volume of the two boxes is not restricted to be constant. The 
particle transfer step to equilibrate the chemical potential in both boxes is carried out 
after each five hundred time steps by comparing their chemical potentials. It is done as 
follows: A particle is chosen randomly from the box where chemical potential is more 
and is removed from it. Correspondingly a particle is introduced into the other with its 
potential energy calculated and the velocity taken from the Boltzmann distribution of 
corresponding temperature. However, care is taken so that the introduced particle is not 
too close to any other particle in the box. With the above three procedures being done 
during the simulation, the two boxes evolve in time in such away that they have same 
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temperature, pressure and the chemical potential after equilibration. Thus the system 
may phase separate into liquid in one box and gas in the other with proper choices of 
initial densities if the temperature of the system is less than the critical temperature. 
The method described above has been tested for the Lennard- Jones fluid initially. The 
phase diagram obtained matched with that of the earlier simulations [H]] validating the 
code we developed. Then the code is used to obtain the LVPDs for the metals described 
above. Time step used for the simulations is 1 femto second and each simulation run 
has typically 5 x 10 5 equilibration and production steps. The chemical potential has 
been evaluated using Widom's test particle insertion method[T7]. In each step 200 test 
particles are inserted and the chemical potential is calculated after each five hundred 
steps. Total number of particles used in the simulation are 1728. The initial densities 
have been chosen so that after equilibration, both the boxes contain a good number of 
particles (few hundreds) so that the averages are reliable and deviations are small. 

3.1. Results and Discussion 

The Temperature (T) Versus Density(p) diagrams for Aluminum, Copper, Sodium and 
Potassium are shown in fig.©. The critical temperature (T c ) and critical density (p c ) 
are obtained by fitting the simulation data using the law of rectilinear diameters 



Critical parameters that we obtained are compared with literature data and 
experiments in Table. (//). In the case of Aluminum, the T c and P c we obtained are higher 
than the literature data. Whereas for copper our data is closer to the experimental 
value than that of Ref. [6 J and Ref.[23]. On the other hand, in the case of Sodium 
and Potassium, our results are significantly improved over those obtained from Morse 
potential[6]. Also, experimental coexistence points are shown in Fig. ([6b) and Fig.([6U) for 
Sodium and Potassium respectively. The figures clearly show a significant improvement 
over the earlier results using Morse potential throughout the phase diagram. However, 
still there is some deviation between our results and experimental coexistence points. 
In the case of Sodium, the coexistence points obtained by Singh et. al.[26] using EAM 
potential are closer to the experimental points at low temperatures. However, the 
asymmetry in experimental phase coexistence curve is not present in that obtained by 
Singh et. al. particularly in the liquid part of the coexistence curve. Whereas the 
asymmetry of the experimental curve is qualitatively maintained by that we obtained 
in this work. 





where pi and p v are liquid and vapor densities. A and B are fitting constants and 
P = 0.33. 



Molecular Dynamic Simulation of Liquid- Vapor Coexistence of Metals Modeled Using Modified Empirical Pa 



Table 2. Critical Point Data 



metal 


T C (K) 


Pc(g/cc) 


P c (GPa) 


Reference 


A 1 

Al 


9643 


0.75 


0.81 


mi • l 

1ms work 




7963 


0.44 


0.35 


Faussurier[T9~] 




8860 


0.28 


0.31 


T ■ l 1 i lr\r\\ 

Likalter|20j 




8387 


0.38 


0.45 


Vinayak[2I] 




8472 


0.79 


0.51 


J.K. Singh[6] 


Uu 






n to 
U./o 


Ims work 




5696 


1.8 


0.11 


All 1 |no 1 

Aleksandrov [23J 




8650 


2.6 


0.95 


J.K. Singh[6] 




7696 


1.93 


0.58 


Experiment |24j 


IN a 


3121 


0.31 


n i fin 

0.109 


rpl • ^ 1 

ims work 




2485 


0.30 


0.025 


Experiment |25j 






n or 
U.ofct 


u.izy 


J.K. omgnjoj 




2462 


0.35 


0.011 


J.K. Singh[26] 


K 


2280 


0.27 


0.037 


This work 




2280 


0.19 


0.016 


Experiment [25] 




3120 


0.28 


0.053 


J.K. Singh[6] 



4. Conclusion 

We proposed a modified form of empirical pair potential for metals and suggested a 
simpler way of obtaining parameters of the potential using the ab-initio cold curve. We 
noted that the cold curve obtained from the modified potential is accurate enough far 
from the equilibrium region for all the metals we worked with in the paper. We obtained 
liquid-vapor phase coexistence curves with the potential using GEMD and noted that 
there is a significant improvement in the phase diagrams over those obtained from Morse 
potential parametrized using equilibrium data. However, we find that still there is some 
deviation from the experimental results in the case of alkali metals. The deviation could 
be artifact of many body effects which cannot be accurately accounted within the pair 
potential formalism. 
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Figure 6. Liquid Vapor Coexistence Curves, (filled squares: Present work), (stars: 
J.K. Singh ct. al. [6]), (Filled circles: J.K. Singh et. al. [26] ), (Crosses: Experiments 
[27]). (a) Aluminum, (b) Copper, (c) Sodium and (d) Potassium. 
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